#' circosPlot for DIABLO
#' 
#' Displays variable correlation among different blocks
#' 
#' \code{circosPlot} function depicts correlations of variables selected with
#' \code{block.splsda} or \code{block.spls} among different blocks, using a
#' generalisation of the method presented in González et al 2012. If
#' \code{ncomp} is specified, then only the variables selected on that component
#' are displayed.
#' 
#' The \code{linkWidth} argument specifies the width of the links drawn.
#' If a vector of length 2 is provided, the smaller value will correspond to
#' a similarity values designated by \code{cutoff} argument, while the 
#' larger value will be used for a link with perfect similarity (1), if any.
#' @param object An object of class inheriting from \code{"block.plsda"},
#' \code{"block.splsda"}, \code{"block.pls"} or \code{"blocks.spls"}.
#' @param comp Numeric vector indicating which component to plot. Default to
#' all
#' @param cutoff Only shows links with a correlation higher than \code{cutoff}
#' @param blocks Character or integer vector indicating which blocks to show.
#' Default to all
#' @param color.Y a character vector of colors to be used for the levels of the
#' outcome
#' @param color.blocks a character vector of colors to be used for the blocks
#' @param color.cor a character vector of two colors. First one is for the
#' negative correlation, second one is for the positive correlation
#' @param var.names Optional parameter. A list of length the number of blocks
#' in \code{object$X}, containing the names of the variables of each block. If
#' \code{NULL}, the colnames of the data matrix are used.
#' @param showIntraLinks if TRUE, shows the correlation higher than the
#' threshold inside each block.
#' @param line if TRUE, shows the overall expression of the selected variables.
#' see examples.
#' @param size.legend size of the legend
#' @param ncol.legend number of columns for the legend
#' @param size.variables size of the variable labels
#' @param size.labels size of the block labels
#' @param legend Logical. Whether the legend should be added. Default is TRUE.
#' @param legend.title String. Name of the legend. Defaults to "Expression".
#' @param linkWidth Numeric. Specifies the range of sizes used for lines linking
#' the correlated variables (see details). Must be of length 2 or 1. Default to c(1). See details.
#' @param ... For object of class \code{block.splsda}, advanced plot parameters:
#' \itemize{
#'  \item \bold{var.adj} Numeric. Adjusts the radial location of variable names in 
#'  units of the arc band width. Positive values push feature names radially 
#'  outward. Default to -0.33. See examples.
#'  \item \bold{block.labels.adj} Numeric between -1 and 1. Adjusts the radial
#'  location of block names radially inward or outward. Default to 0. See examples.
#'  \item \bold{blocks.link} Character vector of blocks. If provided, only correlations
#'  from features of these blocks are shown using links. See examples.
#' }
#' For object of class \code{block.spls}, all listed and advanced arguments
#' passed to the \code{block.splsda} method.
#' @return If saved in an object, the circos plot will output the similarity
#' matrix and the names of the variables displayed on the plot (see
#' \code{attributes(object)}).
#' @author Michael Vacher, Amrit Singh, Florian Rohart, Kim-Anh Lê Cao, Al J Abadi
#' @seealso \code{\link{block.splsda}}, references and
#' http://www.mixOmics.org/mixDIABLO for more details.
#' @references Singh A., Gautier B., Shannon C., Vacher M., Rohart F., Tebbutt
#' S. and Lê Cao K.A. (2016). DIABLO: multi omics integration for biomarker
#' discovery. BioRxiv available here:
#' \url{http://biorxiv.org/content/early/2016/08/03/067611}
#' 
#' mixOmics article:
#' 
#' Rohart F, Gautier B, Singh A, Lê Cao K-A. mixOmics: an R package for 'omics
#' feature selection and multiple data integration. PLoS Comput Biol 13(11):
#' e1005752
#' 
#' González I., Lê Cao K.A., Davis M.J., Déjean S. (2012). Visualising
#' associations between paired 'omics' data sets. \emph{BioData Mining};
#' \bold{5}(1).
#' @keywords regression multivariate
#' @example ./examples/circosPlot-examples.R
#' @importFrom reshape2 dcast
#' @importFrom tidyr gather
#' @importFrom dplyr group_by mutate summarise arrange 
#' @name circosPlot
NULL
#' @export
#' @noRd
circosPlot <- function(object, ...) UseMethod('circosPlot')

#' @keywords Internal
.circosPlot <- 
    function(object,
             comp = 1:min(object$ncomp),
             cutoff,
             color.Y,
             blocks = NULL,
             color.blocks,
             color.cor,
             var.names = NULL,
             showIntraLinks = FALSE,
             line = FALSE,
             size.legend = 0.8,
             ncol.legend = 1,
             size.variables = 0.25,
             size.labels = 1,
             legend = TRUE,
             legend.title = "Expression",
             linkWidth = 1,
             ...)
    {
        
        if (inherits(object, c("block.plsda", "block.splsda"))) {
            indY <- object$indY
            object$variates[indY] <- NULL
            object$loadings[indY] <-  NULL
        }
        
        # to satisfy R CMD check that doesn't recognise x, y and group (in aes)
        Features = Exp = Dataset = Mean = linkColors = chrom = po = NULL
        
        
        figSize = 800
        segmentWidth = 25
        linePlotWidth = 90
        
        ##############################
        ###   networkDiagram_core.R
        ###
        ###   Authors: Michael Vacher (minor changes by Amrit :)
        ###
        ###   Parts of this code has been modified from the original OmicCircos
        ###     package obtained from:
        ###   Ying Hu Chunhua Yan <yanch@mail.nih.gov> (2015). OmicCircos:
        ###     High-quality circular visualization of omics data. v1.6.0.
        ##############################
        
        # check input object
        if (!(inherits(object, c("block.plsda", "block.splsda", "block.pls", "block.spls"))))
            stop("circosPlot is only available for 'block.(s)pls(da)' objects")
        
        if (length(object$X) < 2)
            stop("This function is only available when there are more than 3 blocks
    (2 in object$X + an outcome object$Y)") # so 2 blocks in X + the outcome Y
        if (!is.null(blocks))
        { ## allow selection of blocks to show
            # TODO checks for valid blocks
            
            oblocks <- names(object$X)[-object$indY]
            
            if (is.numeric(blocks))
            {
                blocks <- oblocks[blocks]
            }
            
            drop.blocks <- which(!oblocks %in% blocks)
            if (length(drop.blocks) > 0) {
                object$X <- object$X[-drop.blocks]
                object$variates <- object$variates[-drop.blocks]
                object$loadings <- object$loadings[-drop.blocks]
            }
        }
        
        if (missing(cutoff))
            stop("'cutoff' is missing", call.=FALSE) # so 2 blocks in X + the outcome Y
        
        if(missing(color.Y))
        {
            color.Y = color.mixo(1:nlevels(object$Y))
        } else {
            if(length(color.Y) != nlevels(object$Y))
                stop("'color.Y' must be of length ", nlevels(object$Y))
            
        }
        if (missing(color.blocks))
        {
            color.blocks = brewer.pal(n = 12, name = 'Paired') #why 12?? ANS: bc max allowed n is 12 for this function
            if (length(object$X) > 6) {
                ## if more than 6 X, extend the palette
                ## 2*length(object$X) because we need dark and clear cols for each block
                color.blocks <- colorRampPalette(color.blocks)(2*length(object$X))
            }
            
        } else {
            if(length(color.blocks) != length(object$X))
                stop("'color.blocks' must be of length ", length(object$X))
            
            color.blocks.adj = adjustcolor(color.blocks, alpha.f = 0.5)
            #to get two shades of the same color per block
            
            color.blocks = c(rbind(color.blocks, color.blocks.adj))
            # to put the color next to its shaded color
        }
        
        if(missing(color.cor))
        {
            color.cor = c(colors()[134],  # blue, negative correlation
                          colors()[128])  # pale red, positive correlation
        } else {
            if(length(color.cor) != 2)
                stop("'color.cor' must be of length 2")
        }
        
        ## create unique feature names for circosPlot
        object$X <- mapply(object$X, names(object$X), FUN = function(mat,blockname)
        {
            colnames(mat) <- sprintf("%s_%s", colnames(mat), blockname)
            mat
        }, SIMPLIFY = FALSE)
        loadings.X <- object$loadings
        original.var.names <- lapply(loadings.X, rownames)
        loadings.X <- mapply(loadings.X, names(loadings.X), FUN = function(mat,blockname)
        {
            rownames(mat) <- sprintf("%s_%s", rownames(mat), blockname)
            mat
        }, SIMPLIFY = FALSE)
        object$loadings <- loadings.X
        ## DONE - create unique feature names for circosPlot
        Y = object$Y
        X = object$X
        
        #check var.names
        sample.X = lapply(object$loadings,
                          function(x){1 : nrow(x)})
        if (is.null(var.names))
        {
            var.names.list = original.var.names
        } else if (is.list(var.names)) {
            if (length(var.names) != length(object$loadings))
                stop.message('var.names', sample.X)
            
            if(sum(sapply(1 : length(var.names), function(x){
                length(var.names[[x]]) == length(sample.X[[x]])})) !=
               length(var.names))
                stop.message('var.names', sample.X)
            
            var.names.list = var.names
        } else {
            stop.message('var.names', sample.X)
        }
        
        
        if(any(comp > min(object$ncomp)))
        {
            warning("Limitation to ",min(object$ncomp),
                    " components, as determined by min(object$ncomp)")
            comp[which(comp > min(object$ncomp))] = min(object$ncomp)
        }
        comp = unique(sort(comp))
        
        ## check linkWidth
        invalid.linkWidth <- FALSE
        if (mode(linkWidth) == 'numeric')
        {
            if (length(linkWidth) == 1)
                linkWidth <- rep(linkWidth, 2)
            else if (length(linkWidth) != 2)
                invalid.linkWidth <- TRUE
            
            linkWidth <- sort(linkWidth)
        } else {
            invalid.linkWidth <- TRUE
        }
        
        if (invalid.linkWidth)
            stop("'linkWidth' must be a numeric of length 2 (or 1) specifying ",
                 "the range of widths used for link lines based on similarity measures.", call. = FALSE)
        
        keepA = lapply(object$loadings, function(i)
            apply(abs(i)[, comp, drop = FALSE], 1, sum) > 0)
        cord = mapply(function(x, y, keep){
            cor(x[, keep], y[, comp], use = "pairwise")
        }, x=X, y=object$variates,
        keep = keepA, SIMPLIFY = FALSE)
        
        simMatList = vector("list", length(X))
        for(i in 1:length(cord))
        {
            for(j in 1:length(cord))
            {
                simMatList[[i]][[j]] = cord[[i]] %*% t(cord[[j]])
            }
        }
        simMatList <- lapply(simMatList, function(i) do.call(cbind, i))
        simMat = do.call(rbind, simMatList)
        
        ## merge all similarity measures in one data.frame where rows are samples
        Xdat = as.data.frame(do.call(cbind, X)[, colnames(simMat)])
        ## add Y as another feature
        AvgFeatExp0 <- mutate(.data = Xdat, Y = Y)
        ## create a long data.frame with columns: CELL FEATURE FEATURE_VALUE -Y
        AvgFeatExp0 <- gather(data = AvgFeatExp0, Features, Exp, -Y)
        ## group by Y and feature for averaging and calcuate the mean and SD to 
        ## get a data.frame with columns: Y FEATURES MEAN SD
        AvgFeatExp0 <- group_by(.data = AvgFeatExp0, Y, Features)
        AvgFeatExp0 <- summarise(.data = AvgFeatExp0, Mean = mean(Exp, na.rm = TRUE), SD = sd(Exp,na.rm = TRUE))
        ## add a column as the dataset from which the feature came from:
        AvgFeatExp0$Dataset <- factor(rep(names(X), unlist(lapply(cord, nrow))),
                                      levels = names(X))[match(AvgFeatExp0$Features,colnames(Xdat))]
        # to match Xdat that is reordered in AvgFeatExp0
        featExp <- group_by(.data = AvgFeatExp0 , Dataset, Y)
        featExp <- arrange(.data = featExp , Mean)
        #Generate a circular plot (circos like) from a correlation matrix (pairwise)
        #
        # Args:
        #   simMat: the main correlation matrix.
        #         -> colnames == rownames (pairwise)  values = correlations
        #   featExp: data.frame holding the expression data.
        #   cutoff: minimum value for correlations (<threshold will be ignored)
        #   figSize: figure size
        #   segmentWidth: thickness of the segment (main circle)
        #   linePlotWidth: thickness of the line plot (showing expression data)
        #   showIntraLinks = display links intra segments
        
        
        
        # 1) Generate karyotype data
        chr = genChr(featExp, color.blocks = color.blocks)
        ## ------ re-order blocks based on object$X names
        ## define an ordered factor to re-order blocks so with given X, block orders
        ## are always the same
        Xblocks <- unique(paste0('chr',names(object$X)))
        chr$block <- factor(chr$chrom, levels = Xblocks, ordered = TRUE)
        chr <- chr[order(chr$block),]
        chr$block <- NULL
        
        chr.names = unique(chr$chrom) # paste("chr", 1:seg.num, sep="") 
        # Calculate angles and band positions
        db = segAnglePo(chr, seg=chr.names)
        db = data.frame(db)
        
        # 2) Generate Links
        blocks.link <- list(...)$blocks.link
        links = genLinks(chr, simMat, threshold=cutoff, blocks.link = blocks.link)
        if (nrow(links) < 1)
            warning("Choose a lower correlation threshold to highlight
    links between datasets")
        
        # 3) Plot
        # Calculate parameters
        circleR = (figSize / 2.0) -  segmentWidth - linePlotWidth
        linksR = circleR - segmentWidth
        linePlotR = circleR + segmentWidth
        chrLabelsR = (figSize / 2.0)
        
        # replace chr$name by the ones in var.names (matching)
        # matching var.names.list with object$loadings
        ind.match = match(chr$name, unlist(sapply(object$loadings,rownames)))
        chr$name.user = unlist(var.names.list)[ind.match]
        
        opar1=par("mar")
        par(mar=c(2, 2, 2, 2))
        
        plot(c(1,figSize), c(1,figSize), type="n", axes=FALSE, xlab="",
             ylab="", main="")
        
        #save(list=ls(),file="temp.Rdata")
        # Plot ideogram
        drawIdeogram(R=circleR, cir=db, W=segmentWidth,  show.band.labels=TRUE,
                     show.chr.labels=TRUE, chr.labels.R= chrLabelsR, chrData=chr,
                     size.variables = size.variables, size.labels=size.labels,
                     color.blocks = color.blocks, line = line, ...)
        # Plot links
        if(nrow(links)>0)
            drawLinks(R=linksR, cir=db,   mapping=links,   col=linkColors, lineWidth = linkWidth,
                      drawIntraChr=showIntraLinks, color.cor = color.cor)
        
        # Plot expression values
        cTypes = levels(Y)
        #unique(featExp[,1]) #Get the different disease/cancer types (lines)
        #lineCols = rainbow(nrow(cTypes), alpha=0.5)
        lineCols = color.Y
        #color.mixo(1:nlevels(Y))#color.mixo(match(levels(Y), levels(Y)))
        
        # Fixme: remove this loop and send the whole expr dframe to drawLinePlot
        if(line==TRUE)
        {
            for (i in 1:length(chr.names)){
                seg.name = gsub("chr","",chr.names[i])
                #Get data for each segment
                expr = subset(featExp,featExp$Dataset==seg.name)
                
                expr = dcast(expr, formula = Features ~ Y, value.var="Mean")
                ## changed PAM50 to Y
                expr = merge(expr, chr, by.x="Features", by.y="name")
                expr$po = (as.numeric(expr$chromStart) +
                               as.numeric(expr$chromEnd)) / 2.0
                #' @importFrom dplyr rename
                expr = rename(expr, seg.name = chrom, seg.po = po)
                
                # Reorder columns
                cOrder = c(c(grep("seg.name", colnames(expr)),
                             grep("seg.po", colnames(expr))),c(1:length(cTypes)+1))
                expr = expr[, cOrder]
                
                # Plot data on each sub segment
                subChr = subset(db, db$seg.name == chr.names[i] )
                drawLinePlot(R=linePlotR, cir=subChr,   W=linePlotWidth,
                             lineWidth=1, mapping=expr, col=lineCols, scale=FALSE)
            }
        }
        opar=par("xpd")
        par(xpd=TRUE) # to authorise the legend to be written outside the margin,
        #       otherwise it's too small
        # Plot legend
        if(legend == TRUE)
        {
            # First legeng bottom left corner
            legend(x=5, y = (circleR/4), title="Correlations",
                   c("Positive Correlation", "Negative Correlation"),
                   col = color.cor, pch = 19, cex=size.legend, bty = "n")
            # Second legend bottom righ corner
            if(line==TRUE)
                legend(x=figSize-(circleR/3), y = (circleR/3), title=legend.title,
                       legend=levels(Y),  ## changed PAM50 to Y
                       col = lineCols, pch = 19, cex=size.legend, bty = "n",ncol=ncol.legend)
            # third legend top left corner
            legend(x=figSize-(circleR/2), y = figSize, title="Correlation cut-off",
                   legend=paste("r", cutoff, sep = "="),
                   col = "black", cex=size.legend, bty = "n")
            
            legend(x=-circleR/4, y = figSize, legend=paste("Comp",
                                                           paste(comp,collapse="-")),
                   col = "black", cex=size.legend, bty = "n")
        }
        par(xpd=opar,mar=opar1)# put the previous default parameter for xpd
        
        # use original feature names in similarity matrix if possible
        feat_block <- colnames(simMat)
        # remove anything from last underscore to end of the feature name
        feat <- gsub('(.*)_\\w+', '\\1', feat_block)
        if(!any(duplicated(feat)))
        {
            dimnames(simMat) <- list(feat, feat)
        } else {
            cat("\n")
            cat("adding block name to feature names in the output similarity ")
            cat("matrix as there are similar feature names across blocks.")
            cat("\n")
        }
        return(invisible(simMat))
    }

#' @method circosPlot block.splsda
#' @rdname circosPlot
#' @export
circosPlot.block.splsda <- .circosPlot

#' @method circosPlot block.plsda
#' @rdname circosPlot
#' @export
circosPlot.block.plsda <- .circosPlot

#' @method circosPlot block.spls
#' @rdname circosPlot
#' @param group The grouping factor used when \code{line = TRUE}
#' @param Y.name Character, the name of the \code{Y} block
#' @export
circosPlot.block.spls <- function(object, ..., group = NULL, Y.name = 'Y')
{
    # when a block.spls object is supplied that uses the indY parameter, the name
    # of the object$X component for this dataframe is its proper name
    # for below checks, change it back to "Y"
    if (is.null(object$X$Y)) {
        names(object$X)[object$indY] <- "Y"
    }
    
    if (length(group) != nrow(object$X$Y))
        stop("group must be a factor of length: nrow(object$X$Y) = ", nrow(object$X$Y), "\n")
    object$Y <- factor(group)
    object$keepX <- c(object$keepX, list(Y = object$keepY))
    nX <- length(object$variates)
    
    block.names <- names(object$X)
    block.names[which(block.names == 'Y')] <- Y.name
    
    names(object$X) <-
        names(object$variates) <- 
        names(object$loadings) <- 
        block.names
    .circosPlot(object, ...)
}

#' @method circosPlot block.pls
#' @rdname circosPlot
#' @export
circosPlot.block.pls <- circosPlot.block.spls

## ------------- circosPlot utils
drawIdeogram = function(R, xc=400, yc=400, cir, W,
                        show.band.labels = FALSE,
                        show.chr.labels = FALSE, chr.labels.R = 0,
                        chrData,
                        size.variables,
                        size.labels,
                        color.blocks,
                        line,
                        var.adj = -0.33, ## radial adjustment of feature names
                        block.labels.adj = 0, ## radial adjustment of block label names
                        ...)
{
    # Draw the main circular plot: segments, bands and labels
    chr.po    = cir 
    chr.po[,1]  = gsub("chr","",chr.po[,1]) 
    chr.num     = nrow(chr.po) 
    dat.c     = chrData 
    dat.c[,1] = gsub("chr", "", dat.c[,1]) 
    
    for (chr.i in c(1:chr.num)){
        chr.s  = chr.po[chr.i,1] 
        
        v1 = as.numeric(chr.po[chr.i,2]) 
        v2 = as.numeric(chr.po[chr.i,3]) 
        v3 = as.numeric(chr.po[chr.i,6]) 
        v4 = as.numeric(chr.po[chr.i,7]) 
        
        dat.v = subset(dat.c, dat.c[,1]==chr.s) 
        dat.v = dat.v[order(as.numeric(dat.v[,2])),] 
        for (i in 1:nrow(dat.v)){
            
            #col.v = which(colors()==dat.v[i,5])  #get color index
            #col = colors()[col.v]
            dark.clear = color.blocks#brewer.pal(n = 12, name = 'Paired')
            col.v = which(dark.clear==dat.v[i,5])  #get color index
            col = dark.clear[col.v]
            
            w1 = scale.v(as.numeric(dat.v[i,2]), v1, v2, v3, v4) 
            w2 = scale.v(as.numeric(dat.v[i,3]), v1, v2, v3, v4) 
            
            draw.arc.s(xc, yc, R, w1, w2, col=col, lwd=W) 
            
            if (show.band.labels){
                band.text = as.character(dat.v[i,"name.user"]) 
                
                band.po = ((w1+w2)/2)# - ((w2-w1)/3) #position around the circle
                # print(c(band.po, w1, w2, (w2-w1)/3))
                band.po.in = R + W*var.adj #position on the band (middle)
                draw.text.rt(xc, yc,band.po.in  , band.po , band.text ,
                             cex = size.variables, segmentWidth = W, side="in" )
            }
        } #End for row
        if (show.chr.labels){
            w.m = (v1+v2)/2 
            chr.t = gsub("chr", "", chr.s)
            if(line == TRUE)
            {
                draw.text.rt(xc, yc, chr.labels.R, w.m, chr.t, cex=size.labels,
                             segmentWidth = W, parallel=TRUE)
            } else {
                #put the labels closer to the circle
                draw.text.rt(xc, xc, chr.labels.R, w.m, chr.t, cex=size.labels,
                             segmentWidth = 55 + ceiling(block.labels.adj*55), parallel=TRUE)
            }
        }
    } #End for
}

drawLinks = function(R, xc=400, yc=400, cir, W,
                     mapping=mapping,
                     lineWidth=1, col=rainbow(10, alpha=0.8)[7],  drawIntraChr=FALSE,
                     color.cor = color.cor)
{
    # Draw the links (computed correlation) between features
    chr.po    = cir 
    chr.po[,1]  = gsub("chr","",chr.po[,1]) 
    chr.num     = nrow(chr.po) 
    
    chr.po[,4] = gsub("chr", "", chr.po[,4]) 
    dat.in = mapping 
    dat.in[,1] = gsub("chr", "", dat.in[,1]) 
    dat.in[,4] = gsub("chr", "", dat.in[,4]) 
    
    dat    = dat.in 
    ## line width based on on abs(cor) 
    lineWidth.range = lineWidth
    lineWidths = abs(dat$value)
    ## normlaised so that (cutoff, 1) would map to lineWidth.range
    lineWidths = (lineWidths - min(lineWidths)) / (1 - min(lineWidths))
    lineWidths = lineWidths*(lineWidth.range[2] - lineWidth.range[1]) + lineWidth.range[1]
    for (i in 1:nrow(dat)){
        chr1.s   = dat[i,1] 
        chr2.s   = dat[i,4] 
        po1      = dat[i,2] 
        po2      = dat[i,5] 
        lineWidth = lineWidths[i]
        
        chr1     = which(chr.po[,1]==chr1.s) 
        chr2     = which(chr.po[,1]==chr2.s) 
        
        v1 = as.numeric(chr.po[chr1,2]) 
        v2 = as.numeric(chr.po[chr1,3]) 
        v3 = as.numeric(chr.po[chr1,6]) 
        v4 = as.numeric(chr.po[chr1,7]) 
        
        w1 = scale.v(as.numeric(po1), v1, v2, v3, v4) 
        
        v1 = as.numeric(chr.po[chr2,2]) 
        v2 = as.numeric(chr.po[chr2,3]) 
        v3 = as.numeric(chr.po[chr2,6]) 
        v4 = as.numeric(chr.po[chr2,7]) 
        
        w2 = scale.v(as.numeric(po2), v1, v2, v3, v4) 
        # Set the link width depending on the correlation coefficient
        lwd = abs(as.numeric(dat[i,7])) 
        
        # Set link color
        if (as.numeric(dat[i,7]) < 0.0){
            linkCol = color.cor[2]#colors()[128]  #pale red
        } else {
            linkCol = color.cor[1]#colors()[134]  # blue
        }
        linkCol = add.alpha(linkCol, alpha=0.4) 
        
        if (chr1 == chr2){
            if (drawIntraChr == TRUE){
                draw.link(xc, yc, R, w1, w2, col=linkCol, lwd=lineWidth) 
            }
        } else {
            draw.link(xc, yc, R, w1, w2, col=linkCol, lwd=lineWidth) 
        }
    } ### End for
}

drawLinePlot = function(mapping=mapping, xc=400, yc=400, col.v=3,
                        R, cir,   W, col='black', scale=FALSE, lineWidth=1,
                        background.lines=FALSE,axis.width=1)
{
    # Generate a linear plot around the main ideogram.
    #
    # fixme: the function writes the same line multiple times
    # it needs to be called only once and process the data/segment
    # separately
    chr.po    = cir 
    chr.po[,1]  = gsub("chr","",chr.po[,1]) 
    chr.num     = nrow(chr.po) 
    
    dat.in   = mapping 
    dat.in[,1] = gsub("chr", "", dat.in[,1]) 
    
    # data set for the chromosome
    for (chr.i in 1:chr.num){
        chr.s = chr.po[chr.i,1] 
        chr.s = gsub("chr","",chr.s) 
        dat   = subset(dat.in, dat.in[,1]==chr.s) 
        dat   = dat[order(as.numeric(dat[,2])),] 
        v1 = as.numeric(chr.po[chr.i,2]) ## angle.start
        v2 = as.numeric(chr.po[chr.i,3]) ## angle.end
        v3 = as.numeric(chr.po[chr.i,6]) ## seg.start
        v4 = as.numeric(chr.po[chr.i,7]) ## seg.end
        
        # background line
        if (background.lines){
            draw.arc.pg(xc, yc, v1, v2, R, R+W-5, col=colors()[245]) 
        } else {
            # draw.arc.s(xc, yc, R, v1, v2, col=colors()[245], lwd=axis.width) 
        }
    }
    
    my.R1 = R + W/5 
    my.R2 = R + W - W/5 
    
    ## for the matrix colors
    num.col = ncol(dat[,col.v:ncol(dat)]) 
    num.row = nrow(dat.in) 
    
    if (length(col) == num.col){
        colors = col 
    } else {
        colors = rainbow(num.col, alpha=0.5) 
    }
    
    for (chr.i in 1:chr.num){
        chr.s = chr.po[chr.i,1] 
        chr.s = gsub("chr","",chr.s) 
        
        dat   = subset(dat.in, dat.in[,1]==chr.s) 
        dat   = dat[order(as.numeric(dat[,2])),] 
        #print(head(dat))
        dat.i   = c(col.v:ncol(dat)) 
        dat.m   = dat.in[,dat.i] 
        dat.m   = as.matrix(dat.m) 
        dat.min = min(as.numeric(dat.m), na.rm=TRUE)
        dat.max = max(as.numeric(dat.m), na.rm=TRUE) 
        
        v1 = as.numeric(chr.po[chr.i,2]) 
        v2 = as.numeric(chr.po[chr.i,3]) 
        v3 = as.numeric(chr.po[chr.i,6]) 
        v4 = as.numeric(chr.po[chr.i,7]) 
        
        col.i = 0 
        
        for (j in col.v:ncol(dat)){
            col.i = col.i + 1 
            col   = colors[col.i] 
            
            my.v      = as.numeric(dat[1,j]) 
            dat.i.old = my.v 
            v.old   = scale.v(my.v, my.R1, my.R2, dat.min, dat.max) 
            
            po      = as.numeric(dat[1,2]) 
            w.from  = scale.v(po, v1, v2, v3, v4) 
            
            
            for (k in 1:nrow(dat)){
                
                dat.i = as.numeric(dat[k,j]) 
                
                if (is.na(dat.i)){
                    next 
                }
                
                v    = scale.v(dat.i, my.R1, my.R2, dat.min, dat.max) 
                w.to = scale.v(as.numeric(dat[k,2]), v1, v2, v3, v4) 
                
                if (w.from > 0){
                    draw.line3(xc, yc, w.from, w.to, v.old, v, col=col,
                               lwd=lineWidth)
                }
                
                dat.i.old = dat.i 
                w.from    = w.to 
                v.old     = v 
            } # end the row
        }   # end the col
        if (scale){
            do.scale(xc, yc, dat.min, dat.max, R+W/5, W-2*(W/5)) 
        }
    }     # end the chr/segment
}

genChr =function (expr, bandWidth = 1.0, color.blocks)
{
    # Generate the segments and calculate the
    # unique positions of the bands
    #
    # Args:
    #   expr : dataframe containing the features expression
    # example: colnames(concatFeatExp) "PAM50"    "Features" "Mean"     "SD"
    # "Dataset"
    #   bandWidth: thickness of each band
    #
    # Return:
    #   a data.frame that can be used with segAnglePo
    
    # expr can contains expression data for multiple diseases
    # here, we only use the Chrom and Dataset column and remove duplicates
    keeps = c("Features","Dataset") 
    expr = expr[keeps] 
    expr = unique(expr) 
    chrLengths = data.frame(table(expr$Dataset))  
    rownames(chrLengths) = chrLengths[,1] 
    chrLengths[,1] = NULL 
    colnames(chrLengths) = c( "Freq") 
    chrLengths[, "Count"] = rep(0, nrow(chrLengths)) 
    
    # Last column contains the bands' color
    # Create a color scheme
    #dark = c("brown3","darkgoldenrod","antiquewhite3","steelblue3")
    #clear = c("brown1","darkgoldenrod1","antiquewhite1","steelblue1")
    dark.clear = color.blocks#brewer.pal(n = 12, name = 'Paired')
    dark = dark.clear[seq(2, length(dark.clear), by = 2)]
    clear = dark.clear[seq(1, length(dark.clear), by = 2)]
    chrColScheme = data.frame(dark, clear)
    n_datasets = length(unique(expr$Dataset))
    chrColScheme = chrColScheme[c(1:n_datasets),]
    rownames(chrColScheme) = levels(factor(expr$Dataset))#alphabetical order
    
    seg.out = c() 
    for (i in 1:nrow(expr)){
        chrName    = paste("chr", as.character(expr[i,'Dataset'][[1]]), sep="") 
        dType = as.character(expr[i,'Dataset'][[1]]) 
        pStart = chrLengths[dType,'Count'] * bandWidth 
        pStop = chrLengths[dType,'Count'] * bandWidth + bandWidth  
        chrLengths[dType,'Count'] = chrLengths[dType,'Count'] + 1 
        fName = as.character(as.matrix(expr[i,'Features']))
        # added as.character() by amrit
        # Assign colors
        if (chrLengths[dType,'Count'] %% 2 == 0){
            chrCol = chrColScheme[dType,]$clear 
        } else{
            chrCol = chrColScheme[dType,]$dark 
        }
        seg.out = rbind(seg.out, c(chrName, pStart, pStop,  fName, chrCol)) 
    }
    
    # Use the same names than in omicCircos
    colnames(seg.out) = c("chrom", "chromStart", "chromEnd", "name", "color") 
    seg.out = as.data.frame(seg.out) 
    
    return(seg.out) 
}

#' @importFrom reshape2 melt 
genLinks = function(chr, simMat, threshold, blocks.link = NULL)
{
    
    # to satisfy R CMD check that doesn't recognise x, y and group (in aes)
    Var1=Var2=chrom=NULL
    
    
    # Generates the links corresponding to pairwise correlations
    #
    # Args:
    #   chr: ideogram structure (generated from genChr)
    #   simMat: main correlation matrix
    #   threshold: minimum correlation value
    #
    # Return:
    #   the links data (see omicsCircos doc)
    #
    linkList = c() 
    # Remove matrix diagonal and the the mat in a list
    linkList = subset(melt(simMat), Var1!=Var2 ) 
    # Remove links below the threshold
    linkList = subset(linkList, abs(linkList$value) >= threshold) 
    
    ## Al: if we want links only to certain block
    if (!is.null(blocks.link)) {
        only.chr <- paste0('chr', blocks.link)
        only.names <- chr[chr$chrom == only.chr,]$name
        linkList <- linkList[linkList$Var2 %in% only.names,]
    }
    
    #First merge
    #' @importFrom dplyr rename
    linkList = rename(linkList, feat1=Var1, feat2=Var2)
    # CHANGED BY AMRIT
    linkList = merge(linkList, chr, by.x="feat1", by.y="name") 
    # Set the position in the middle of the band
    linkList$po1 = (as.numeric(linkList$chromStart)
                    + as.numeric(linkList$chromEnd)) / 2.0
    #' @importFrom dplyr rename
    linkList = rename(linkList, chr1=chrom)   # CHANGED BY AMRIT
    keeps = c("feat1","feat2","value","chr1","po1") 
    linkList = linkList[keeps] 
    
    #Second merge
    linkList = merge(linkList, chr, by.x="feat2", by.y="name") 
    linkList$po2 = (as.numeric(linkList$chromStart)
                    + as.numeric(linkList$chromEnd)) / 2.0
    #' @importFrom dplyr rename
    linkList = rename(linkList, chr2=chrom)   # CHANGED BY AMRIT
    keeps = c("chr1","po1","feat1","chr2","po2","feat2","value") 
    linkList = linkList[keeps] 
    
    return(linkList) 
}

bezierCurve = function(x, y, n=10)  {  
    outx = NULL  
    outy = NULL   
    i = 1	
    for (t in seq(0, 1, length.out=n))		{		
        b = bez(x, y, t)		
        outx[i] = b$x		
        outy[i] = b$y 		
        i = i+1		
    } 	
    return (list(x=outx, y=outy))	
} 

##
bez = function(x, y, t)	{	
    outx = 0	
    outy = 0	
    n = length(x)-1	
    for (i in 0:n)		{		
        outx = outx + choose(n, i)*((1-t)^(n-i))*t^i*x[i+1]		
        outy = outy + choose(n, i)*((1-t)^(n-i))*t^i*y[i+1]		
    } 	
    return (list(x=outx, y=outy))	
}

###########################################
# one value : from a to b
scale.v = function(v, a, b, min.v, max.v) {
    v = v-min.v 
    v = v/(max.v-min.v) 
    v = v*(b-a) 
    v+a
}

### draw.link
draw.link = function(xc, yc, r, w1, w2, col=col, lwd=lwd) {
    # for translocation
    w3  = (w1+w2)/2 
    w1  = w1/360*2*pi 
    w2  = w2/360*2*pi 
    w3  = w3/360*2*pi 
    x0  = xc+r*cos(w1) 
    y0  = yc-r*sin(w1) 
    x1  = xc+r*cos(w2) 
    y1  = yc-r*sin(w2) 
    x = c(x0,xc,xc,x1) 
    y = c(y0,yc,yc,y1) 
    points(bezierCurve(x,y,60), type="l", col=col, lwd=lwd, lend="butt")
}

### draw.link2
draw.link2 = function(xc, yc, r, w1, w2, col=col, lwd=lwd) {
    # for translocation
    w3  = (w1+w2)/2 
    w1  = w1/360*2*pi 
    w2  = w2/360*2*pi 
    w3  = w3/360*2*pi 
    x2  = xc+r/2*cos(w3) 
    y2  = yc-r/2*sin(w3) 
    x0  = xc+r*cos(w1) 
    y0  = yc-r*sin(w1) 
    x1  = xc+r*cos(w2) 
    y1  = yc-r*sin(w2) 
    x = c(x0, x2, x2, x1) 
    y = c(y0, y2, y2, y1) 
    points(bezierCurve(x,y,60), type="l", col=col, lwd=lwd, lend="butt")
}
###

###
draw.link.pg = function(xc, yc, r, w1.1, w1.2, w2.1, w2.2, col=col, lwd=lwd) {
    w1 = w1.1 
    w2 = w2.2 
    w3  = (w1+w2)/2 
    w1  = w1/360*2*pi 
    w2  = w2/360*2*pi 
    w3  = w3/360*2*pi 
    x0  = xc+r*cos(w1) 
    y0  = yc-r*sin(w1) 
    x1  = xc+r*cos(w2) 
    y1  = yc-r*sin(w2) 
    x = c(x0,xc,xc,x1) 
    y = c(y0,yc,yc,y1) 
    bc1 = bezierCurve(x,y,60) 
    
    ang.d = abs(w1.1-w1.2) 
    pix.n = ang.d * 10 
    if (pix.n < 10){
        pix.n = 10 
    }
    
    ang.seq = rev(seq(w1.1,w1.2,length.out=pix.n)) 
    ang.seq = ang.seq/360*2*pi 
    
    fan.1.x = xc + cos(ang.seq) * r 
    fan.1.y = yc - sin(ang.seq) * r 
    
    ######################################################
    w1 = w1.2 
    w2 = w2.1 
    w3  = (w1+w2)/2 
    w1  = w1/360*2*pi 
    w2  = w2/360*2*pi 
    w3  = w3/360*2*pi 
    x0  = xc+r*cos(w1) 
    y0  = yc-r*sin(w1) 
    x1  = xc+r*cos(w2) 
    y1  = yc-r*sin(w2) 
    x = c(x0,xc,xc,x1) 
    y = c(y0,yc,yc,y1) 
    bc2 = bezierCurve(x,y,60) 
    
    ang.d = abs(w2.1-w2.2) 
    pix.n = ang.d * 10 
    if (pix.n < 10){
        pix.n = 10 
    }
    
    ang.seq = rev(seq(w2.1,w2.2,length.out=pix.n)) 
    ang.seq = ang.seq/360*2*pi 
    
    fan.2.x = xc + cos(ang.seq) * r 
    fan.2.y = yc - sin(ang.seq) * r 
    
    polygon(c(bc1$x, fan.2.x, rev(bc2$x), rev(fan.1.x)),
            c(bc1$y, fan.2.y, rev(bc2$y), rev(fan.1.y)),
            fillOddEven=TRUE, border=col, col=col, lwd=lwd) 
}

###
draw.point.w = function(xc, yc, r, w, col=col, cex=cex){
    w = w/360*2*pi 
    x = xc+r*cos(w) 
    y = yc-r*sin(w) 
    points(x, y, pch=20, col=col, cex=cex) 
}

###
draw.text.w = function(xc, yc, r, w, n, col="black", cex=1){
    w = w%%360 
    w = w/360*2*pi 
    x = xc+r*cos(w) 
    y = yc-r*sin(w) 
    text(x,y,labels=n, col=col, cex=cex) 
}

###
draw.text.rt = function(xc, yc, r, w, n, col="black", cex=1, side="out",
                        segmentWidth=20, parallel=FALSE){
    w     = w%%360 
    the.o = w 
    
    the.w = 360-w 
    w     = w/360*2*pi 
    x     = xc+r*cos(w) 
    y     = yc-r*sin(w) 
    
    
    num2  = (segmentWidth*2)/2.0 
    b = the.w
    if (side=="out"){
        if (the.w <= 90 ){
            the.pos = 4 
            if (parallel == TRUE){
                the.w = the.w -90 # 180 
            }
        } else if (the.w > 90 & the.w <= 180) {
            if (parallel == TRUE){
                the.w = the.w -90 # 180 
            }
            else {
                the.w = the.w + 180 
            }
            the.pos = 2 
        } else if (the.w > 180 & the.w <= 270){
            the.w = the.w%%180 
            if (parallel == TRUE){
                the.w = the.w -90 # 180 
            }
            the.pos = 2 
        } else if (the.w > 270 & the.w <= 360){
            the.pos = 4 
            if (parallel == TRUE){
                the.w = the.w + 90 
            }
        }
        
        if (the.pos==2){
            x = x+num2 
        }
        if (the.pos==4){
            x = x-num2 
        }
    }
    
    if (side=="in"){
        if (the.w <= 90 ){
            the.pos = 4 
        } else if (the.w > 90 & the.w <= 180) {
            the.w = the.w + 180 
            the.pos = 2 
        } else if (the.w > 180 & the.w <= 270){
            the.w = the.w%%180 
            the.pos = 2 
        } else if (the.w > 270 & the.w <= 360){
            the.pos = 4 
        }
        
        if (the.pos==2){
            x = x+segmentWidth 
        }
        if (the.pos==4){
            x = x-segmentWidth 
        }
    }
    
    
    text(x, y, adj=0, offset=1, labels=n, srt=the.w,
         pos=the.pos, col=col, cex=cex) 
}

###strokeLine2
draw.line = function (xc, yc, w, l1, l2, col=col, lwd=lwd, lend=1) {
    w  = (w/360)*2*pi 
    x1 = xc+l1*cos(w) 
    y1 = yc-l1*sin(w) 
    x2 = xc+l2*cos(w) 
    y2 = yc-l2*sin(w) 
    segments(x1, y1, x2, y2, col=col, lwd=lwd, lend=lend) 
}

###strokeLine3
draw.line2 = function (xc, yc, w, r, l, col=col, lwd=lwd){
    line_w   = l 
    theangle = w 
    l1       = r 
    theangle = (theangle/360)*2*pi 
    x0       = xc+l1*cos(theangle) 
    y0       = yc+l1*sin(theangle) 
    w1       = 45/360*2*pi 
    x1 = xc + sin(w1) * (x0) 
    y1 = yc + cos(w1) * (y0) 
    x2 = xc - sin(w1) * (x0) 
    y2 = yc - cos(w1) * (y0) 
    segments(x1, y1, x2, y2, col=col, lwd=lwd, lend="butt") 
}

###strokeLine by two angles
draw.line3 = function (xc, yc, w1, w2, r1, r2, col=col, lwd=lwd){
    theangle1 = w1 
    theangle2 = w2 
    l1        = r1 
    l2        = r2 
    
    theangle1 = (theangle1/360)*2*pi 
    x1        = xc+l1*cos(theangle1) 
    y1        = yc-l1*sin(theangle1) 
    
    theangle2 = (theangle2/360)*2*pi 
    x2        = xc+l2*cos(theangle2) 
    y2        = yc-l2*sin(theangle2) 
    
    segments(x1, y1, x2, y2, col=col, lwd=lwd, lend="butt") 
}

### plot fan or sector that likes a piece of doughnut (plotFan)
draw.arc.pg = function (xc, yc,
                        w1, w2, r1, r2, col="lightblue", border="lightblue", lwd=0.01
){
    
    ang.d = abs(w1-w2) 
    pix.n = ang.d * 10 
    if (pix.n < 10){
        pix.n = 10 
    }
    
    ang.seq = rev(seq(w1,w2,length.out=pix.n)) 
    ang.seq = ang.seq/360*2*pi 
    
    fan.i.x = xc + cos(ang.seq) * r1 
    fan.i.y = yc - sin(ang.seq) * r1 
    
    
    fan.o.x = xc + cos(ang.seq) * r2 
    fan.o.y = yc - sin(ang.seq) * r2 
    
    polygon(c(rev(fan.i.x), fan.o.x ), c(rev(fan.i.y), fan.o.y),
            fillOddEven=TRUE, border=border, col=col, lwd=lwd, lend=1)
    
}

draw.arc.s = function (xc, yc, r, w1, w2, col="lightblue", lwd=1, lend=1){
    # Draw circular arcs for the main ideogram)
    # s = simple
    # r = radius
    ang.d = abs(w1-w2) 
    pix.n = ang.d * 5 
    if (pix.n < 2){
        pix.n = 2 
    }
    
    ang.seq = rev(seq(w1,w2,length.out=pix.n)) 
    ang.seq = ang.seq/360*2*pi 
    
    fan.i.x = xc + cos(ang.seq) * r 
    fan.i.y = yc - sin(ang.seq) * r 
    ## lend=0(round)  #lend=1(butt)  lend=2(square)
    lines(fan.i.x, fan.i.y, col=col, lwd=lwd, type="l", lend=lend) 
    #points(fan.i.x, fan.i.y, col=col, lwd=lwd, type="l", lend=lend) 
}

#########################################
## segment to angle and position
## segAnglePo
#########################################

# get angle if given seg and position
# seg should be ordered by user
segAnglePo = function (seg.dat=seg.dat, seg=seg, angle.start=angle.start,
                       angle.end=angle.end){
    
    if (missing(angle.start)){
        angle.start = 0 
    }
    if (missing(angle.end)){
        angle.end = 360 
    }
    ## check data.frame?
    colnames(seg.dat) = c("seg.name","seg.Start","seg.End","name","gieStain") 
    
    ## get length of the segomosomes
    seg.l   = c() 
    seg.min = c() 
    seg.sum = c() 
    seg.s   = 0 
    seg.num   = length(seg) 
    seg.names = seg 
    
    ########################################################
    ########################################################
    for (i in 1:seg.num){
        seg.n = seg.names[[i]] 
        
        dat.m = subset(seg.dat, seg.dat[,1]==seg.n) 
        seg.full.l = max(as.numeric(dat.m[,"seg.End"])) 
        seg.full.m = min(as.numeric(dat.m[,"seg.Start"])) 
        seg.l      = cbind(seg.l, seg.full.l) 
        seg.min    = cbind(seg.min, seg.full.m) 
        seg.s      = seg.s + seg.full.l 
        seg.sum    = cbind(seg.sum, seg.s) 
    }
    
    ## initial parameters
    gap.angle.size = 2 
    seg.angle.from = angle.start + 270 
    
    seg.full.l  = sum(as.numeric(seg.l)) 
    angle.range = angle.end - angle.start 
    cir.angle.r = (angle.range - seg.num * gap.angle.size)/seg.full.l 
    
    out.s     = c() 
    l.old     = 0 
    gap.angle = 0 
    for (i in 1:seg.num){
        seg.n = seg.names[[i]] 
        dat.m = subset(seg.dat, seg.dat[,1]==seg.n) 
        len   = seg.sum[i] 
        w1    = cir.angle.r*l.old + gap.angle 
        w2    = cir.angle.r*len   + gap.angle 
        out.s     = rbind(out.s, c(seg.n, w1+seg.angle.from, w2+seg.angle.from,
                                   l.old, len, seg.min[i], seg.l[i]))
        gap.angle = gap.angle + gap.angle.size 
        l.old     = len 
    }
    
    colnames(out.s) = c("seg.name","angle.start", "angle.end", "seg.sum.start",
                        "seg.sum.end","seg.start", "seg.end")
    return(out.s) 
}


### start do.scale
do.scale = function(xc=xc, yc=yc, dat.min=dat.min, dat.max=dat.max,
                    R=R, W=W, s.n=1, col="blue"){
    dat.m   = round((dat.min+dat.max)/2, s.n) 
    dat.min = round(dat.min, s.n) 
    dat.max = round(dat.max, s.n) 
    y1      = yc + R  
    y2      = yc + R + W/2 
    y3      = yc + R + W 
    x1      = xc - W/20 
    x2      = x1 - (W/20)*1.2 
    x3      = x1 - (W/20)*3 
    segments(x1, y1, x1, y3, lwd=0.01, col=col) 
    segments(x1, y1, x2, y1, lwd=0.01, col=col) 
    segments(x1, y2, x2, y2, lwd=0.01, col=col) 
    segments(x1, y3, x2, y3, lwd=0.01, col=col) 
    text(x3, y1, dat.min, cex=0.2, col=col) 
    text(x3, y2, dat.m,   cex=0.2, col=col) 
    text(x3, y3, dat.max, cex=0.2, col=col)
}
### end do.scale

## Add an alpha value to a colour
add.alpha = function(col, alpha=1){
    if(missing(col))
        stop("Please provide a vector of colours.")
    apply(sapply(col, col2rgb)/255, 2, 
          function(x) 
              rgb(x[1], x[2], x[3], alpha=alpha))  
}
###-------------------------------------------------------------------------
